/*! \file Example how to simulate a NaCl solution
 *  \example widom-example.C
 *
 * This example program calculates the excess chemical
 * potential of NaCl in an aqueous solution using Widom's
 * particle insertion method. An output POVRAY file will
 * be generated. Note: No equilibration run is incorporated.
 */
#include <iostream>
#include "../io.h"
#include "../analysis.h"
#include "../container.h"
#include "../potentials.h"
typedef pot_test T_pairpot;
#include "../markovmove.C"
#include "../analysis.h"

using namespace std;

int main() {
  cell cell(70.);                       // We want a spherical cell
  canonical nvt;                        // Use the canonical ensemble
  pot_setup cfg;                        // Setup pair potential - default values
  interaction<T_pairpot> pot(cfg);      // Functions for interactions
  saltmove sm(nvt, cell, pot);          // Class for salt movements
  widom<T_pairpot> widom(cell, pot,
      particle::NA, particle::CL);      // Class for Widom particle insertion

  group salt;                                   // Group for mobile ions
  salt.add( cell, particle::NA, 10 );           // Insert sodium ions
  salt.add( cell, particle::CL, 10 );           // Insert chloride ions
  systemenergy sys(pot.energy(cell.p));         // Track system energy

  for (int macro=0; macro<10; macro++) {        // Markov chain
    for (int micro=0; micro<1e2; micro++) {
      sm.move(salt);                            // Displace salt particles
      sm.adjust_dp(40,50);                      // Stride to 40-50% acceptance
      sys+=sm.du;                               // Sum system energy changes
      widom.insert(10);                         // Widom particle insertion analysis
    }
    sys.update(pot.energy(cell.p));             // Update system energy
  }
  cout << cell.info() << sys.info() << sm.info() << widom.info();
}

